Density Functional Theory of Inhomogeneous Liquids: III. Liquid- Vapor Nucleation 
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The process of nucleation of vapor bubbles from a superheated liquid and of liquid droplets 
from a supersaturated vapor is investigated using the Modified-Core van der Waals model Density 
Functional Theory (Lutsko, JCP 128, 184711 (2008)). A novel approach is developed whereby 
nucleation is viewed as a transition from a metastable state to a stable state via the minimum 
free energy path which is identified using the nudged elastic-band method for exploring free energy 
surfaces. This allows for the unbiased calculation of the properties of sub- and super-critical clusters, 
as well as of the critical cluster. For Lennard- Jones fluids, the results compare well to simulation 
and support the view that even at high supersaturation nucleation proceeds smoothly rather than 
via spinodal-like instabilities as has been suggested recently. 
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I. INTRODUCTION 

The liquid-vapor phase transition is of fundamental 
importance in many areas of physics, chemistry and 
engineering. For this reason, the process of homoge- 
neous nucleation of bubbles in a superheated fluid and 
of droplets in a supersaturated vapor has a long history 
going all the way back to the fundamental paper of van 
der Waals [EH!!. 

Recent years have seen tremendous 
progress in the development of simulation methods that 
circumvent the time and size limitations of simulations to 
allow for the study of rare events, such as homogeneous 
nucleation 

UBS 01- In parallel, various approaches to 
a theoretical description of nucleation have been formu- 
lated with varying degrees of success]!, H, EH- However, 
the theoretical description of nucleation is difficult be- 
cause it is fundamentally a nonequilibrium process: the 
initial and final states are metastable and stable, respec- 
tively, but the transition between them involves a se- 
quence of unstable configurations making the use of equi- 
librium statistical mechanics difficult. 

Density Functional Theory (DFT) is the modern real- 
ization of van der Waals' approach to inhomogeneous flu- 
ids. The fundamental quantity in DFT is the ensemble- 
averaged local density, p(r) in the presence of an ex- 
ternal field, <j)(r) at fixed temperature, T and chemical 
potential, p. It can be proven that there is a one to 
one relationship between applied fields and the result- 
ing density profiles. Furthermore, it can also be shown 
that there exists a functional, F[n] such that the quan- 
tity fl[n] = F[n) + J n(r) (<f>(r) — p) dr is minimized when 
n(r) is the equilibrium density profile, p(r) correspond- 
ing to the applied field [TTI. Then, tl[p] is the grand 
potential for the given field, temperature and chemical 
potential. Since the theory only gives meaning to den- 
sity functions which minimize f2[/i], the question is how 
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to extract information from this formalism in the case of 
nucleation in the absence of an external field (or, perhaps, 
at constant gravitational field). 

Intuitively, one imagines that if the functional Q[p] 
gives the free energy for any density profile p(r) in the 
presence of the appropriate stabilizing field, then, in the 
absence of such a field, it should give information about 
the energetic cost of moving from a stable state to the 
state p(r). For closed systems, this intuition has been 
formalized in recent years in the development known 
as dynamical density functional theory [3 EH EH and, 
more generally, it underlies phase-field models for both 
open and closed systems (see, e.g. ref.[l6| for a recent 
example). For all of these models, when fluctuations 
are included, detailed balance implies that the proba- 
bility of observing a density profile p(r) is proportional 
to e~^ L \-P\l a where L[p] = F[p] for closed systems and 
L[p] = Q[p] for open systems and where a characterizes 
the strength of fluctuations. One then expects that an 
approximation to the most likely path from one state to 
another in an open system will be the shortest path in- 
volving the smallest free energy barriers which is to say 
the Minimal Free Energy Path (MFEP) as defined by 
fi[p]. This of course underlies the classical theory of nu- 
cleation and is the view adopted in the present work. 

The simplest method to approximate the MFEP 
is to define a class of parametrized profiles, p(r) — 
p(r; R,w, ...) where R is the radius of the bubble or 
droplet, w is its width and the notation indicates that 
there may be additional parameters as well. Then, the 
functional F[p] becomes simply a function of the param- 
eters R, w, ... and one can define a path by minimizing 
it for each fixed value of R with respect to the remaining 
parameters. The problem, of course, is that one restricts 
the space of possible density profiles and may exclude the 
most likely profile. A second possibility is to introduce 
some auxiliary constraint that defines what is meant by 
a bubble (or droplet) of a given size. For example, such 
a constraint for a bubble could be J (pi — /»(r)) dr = AN 
which defines a profile to be of size AN if the total num- 
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ber of atoms is AN less than in the uniform fluid. For 
reasons described elsewhere [171], this particular example 
of a constraint is not really useful but more realistic con- 
straints have been constructed [§, [Tol. IT§]. Note however 
the philosophy of this approach: a path through density- 
function space is defined by the parameter A7V (or what- 
ever parameters enter the constraint). This is not so 
different in spirit from the use of parametrized profiles 
as both methods effectively reduce the dimensionality of 
the problem. However, in the case of gas bubbles, it 
has been shown using both numerical calculations and 
analytic models that the constraint approach can fail to 
yield a stable profile for a fixed value of A7V (or what- 
ever parameters are used) even though the underlying 
free energy landscape is well behaved[l7(. Hence, while 
it is certainly more general than the use of parametrized 
profiles, the constraint method is not a robust approach 
to the determination of the MFEP. 



In fact, the search for minimum energy paths over some 
energy landscape is a problem that occurs in many ap- 
plications such as the determination of chemical reaction 
paths from ab initio calculations and the determination 
of transitions between cluster structures [lj|. In recent 
years, several techniques have been developed for deter- 
mining the MFEP. Here, one such method - the Nudged 
Elastic Band (NEB) [20I l2l| - is applied to the problem of 
nucleation. This allows for a completely unbiased and ro- 
bust determination of the MFEP between two metastable 
states. 



In this paper, the details of droplet and bubble nucle- 
ation are calculated for the Lennard- Jones fluid. Because 
the structure of the liquid-vapor interface is a balance 
between bulk free energy and surface tension, it is im- 
portant that a useful theory be able to describe both 
quantities with quantitative accuracy. The calculations 
presented here therefore use the Modified-Core van der 
Waals (MC-VDW) model DFT[H which gives a quan- 
titatively accurate description of the planar liquid-vapor 
interface, of the structure of the LJ fluid near a hard 
wall and of the L J fluid confined to a slit pore 22j . In 



the next section, the model is described and the details 
of the application of the NEB method to the problem 
of bubble nucleation is outlined. The results of the cal- 
culations are presented in Section III. There evidence is 
given of the robustness of the method and direct, quan- 
titative comparison is made to the results of computer 
simulation. The MFEP for nucleation of droplets from a 
supersaturated vapor and of bubbles from a supercooled 
liquid are described. The last Section gives a summary 
of the results obtained and their implication concerning 
recent claims of non-standard pathways in liquid-vapor 
nucleation. 



II. THEORY 
A. Density Functional Theory 

The model free energy functional used in my calcula- 
tions is the Modified-Core van der Waals model[22] which 
is written as a sum of three contributions, 

F[p] = F id [p] + F hs \p] + F core [p] + Ft ai i \p\. (1) 

The first contribution is the ideal gas term which is given 
by 



Fu[p] 



00) log 00)) - p(r)) dr. 



(2) 



Next is a hard-sphere contribution, Fh s [p], for which 
the "White Bear" Fundamental Measure Theory (FMT) 
model was used [23l , 1241 a long with the Barker-Hendersen 
hard-sphere diameter[l2. 25]. The third contribution, the 
"core correction" F core [p] , is similar to a FMT model but 
is constructed so that the total free energy functional re- 
produces a given equation of state in the bulk phase as 
well as certain other conditions concerning the direct cor- 
relation function in the bulk fluid 22]. The final term is 
a mean-field treatment of the long-range attraction, 



F t aii[p}= J ©0i2 - d)p(Y 1 )p(Y 2 )v(r 12 )d 



ridr 2 



(3) 



where Q(x) is the step function, d is the Barker- 
Henderson hard-sphere diameter and v(r) is the pair po- 
tential. In most of the calculations described below, the 
potential is a truncated and shifted Lennard- Jones (LJ) 
interaction, 



v(r) 



VLj{r) 



VLj{r c ) 




r < r c 
r > r c 



(4) 



where VLj(r) = 4e ^(~) — (~) J is the untruncated LJ 

potential. (In one case, comparison is made to simula- 
tions performed with an unshifted potential and in that 
case the calculations were performed with a truncated, 
but unshifted potential.) The DFT model requires as in- 
put the bulk equation of state. Since the object of the 
calculations was to model the LJ system as accurately 
as possible, the empirical equation of state of Johnson, 
Zollweg and Gubbins [26| was used. 

B. Determining the MFEP 

The NEB method is a chain-of-states description of 
the MFEP. A path in density space is described by a 
collection of profiles, {p a (r)}al - To be concrete, con- 
sider the problem of the nucleation of bubbles in a su- 
perheated liquid. Then, the initial state is the uniform 
liquid, //°) (r) = pi where p\ is the bulk liquid density 
determined by the temperature and chemical potential. 
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The subsequent images are initialized to a guess of the 
MFEP: for example, a sequence of hyperbolic tangents 
with increasing radii. Of course, if one simply tries to 
minimize the total energy of the path, X^l/ ^[p^- 1 ] then 
nothing is learned since the images will eventually con- 
verge to one of the two attractors, the uniform liquid or 
the uniform vapor, depending on whether their initial val- 
ues are smaller than, or larger than, the critical cluster. 
The starting point of the NEB method is the addition 
of fictitious elastic forces between the images in order 
to force them to remain evenly spaced along the MFEP. 
The following discussion is divided into two parts: first 
the general implementation of the NEB is described after 
which the specialization to spherical symmetry is given. 

The basic element needed to implement the method 
is a definition of a scalar product in density space. The 
natural choice is, for two real functions /(r) and <?(r), 



f*9= / f(r)g(r)de 



(5) 



Note that the limits of the integral are not specified: for 
most purposes here, the product will only be evaluated 
for functions having finite support so that some large 
volume will suffice. Using this definition, the distance 
between two profiles, p^'(r) and p^\v) is defined as 



(6) 



which is the usual L 2 -norm. The idea behind the NEB is 
to minimize the free energy functionals, £l[p l ] , in the man- 
ifold orthogonal to the current estimate of the MFEP and 
to move them along the estimated MFEP using fictitious 
elastic forces to maintain an even spacing of the images. 
To this end, the critical element is the estimation of the 
tangent to the MFEP at each density image for which 



the algorithm of ref. [2fJ was used. This involves the 
neighboring images and their free energies. For example, 
if Sllp^- 1 }] < ft[pW] < Cl[p( i+ % then the tangent to the 



image p 



W called i( l 



is 



t«(r) =p (l+1) (r) -p (i) (r) 



(7) 



and the normalized tangent, £W(r) = i (i) (r)/ (*W * t&). 
If the inequalities are reversed, the tangent is in the di- 



rection p' 



,(*-!) 



For non-monotonic neighbors, the 



heuristic is given in ref. [20| . The NEB method then con- 
sists of finding a configuration that gives zero NEB-force. 
Let the "force" due to the actual free-energy surface be 

solving 



Then the NEB method consists of 



= f 1(!) (r) + #»(r) ( 



w 



p 



(i-1) 



p v ' - p 

(8) 

where F^(r) = jrW ( r ) - £« ( r ) (t» * JF«) is the com- 
ponent of the thermodynamic force orthogonal to the 
tangent vector and k is the spring constant. 

The specialization to spherical geometry is made by 
noting that if the density is spherically symmetric, p(r) = 
p(r), then the corresponding thermodynamic forces, J-(r) 



and J~s{ r ) 
chain rule, 



dp(r) 



are related by the functional 



d/3n[ P ^} dpi?) 

dp(r') dp(r) 
d(3VL[ P ^ 



dr' 



(9) 
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It therefore follows that 



tW( r ) =p { - i+1 \r) -p {{) {r) 



= I t {i) {r)T {i) {v)dv 
t^\r)^\r)dr 



r 



(10) 



p (0 _ p (i-D 



and so forth. Note that this differs somewhat from the 
more heuristic scheme described in ref. 17] . The present 
approach, starting from the general case and specializ- 
ing to spherical symmetry, is more systematic and yields 
lower free energies than did the earlier implementation. 



A final refinement is the use of a so-called "climbing 
image" [21] . This is an image, say image P ^ for which the 
sign of the component of T^ c ' along the path is reversed, 
instead of eliminating this component, thus causing it to 
climb towards the local maximum. For this image, no 
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spring forces are applied so that its behavior is governed 
by 



-Tl 



(r) + (4 c) (r)-^ (c) W) 



= 0. 



(11) 



(Note, however, that the images on either side of the 
climbing image still feel spring forces connecting them 
to the climbing image: there is no analogy of Newton's 
third law in this case.) The climbing image proves very 
effective in determining the saddle point, which is to say 
in this problem, the critical cluster. 

So far, nothing definite has been said about the end- 
point of the chain of states. The initial point is the uni- 
form liquid, but the end cannot be taken to be the uni- 
form gas because the distance between the uniform gas 
and any finite bubble is infinite. In an initial implemen- 
tation of this method [I?) , the end point was taken to be 
a sigmoidal profile with fixed, large radius and with the 
width adjusted to minimize the free energy. The idea was 
that the profiles of large bubbles would be basically sig- 
moidal and any error introduced by the approximation 
would be insignificant. However, this is somewhat unsat- 
isfying as it potentially biases the chain of states away 
from the MFEP. In the present work, a more elegant ap- 
proach was taken. Note that in a part of the chain of 
states where the free energy is monotonically decreasing, 
as it is for clusters larger than the critical cluster, the 
tangent vector for the image i is determined by the im- 
ages i — 1 and i. The only role played by image i + 1 is in 
the calculation of the spring force on image i. However, if 
this image were excluded completely, the spring between 
i and i— 1 would cause these images to converge. To avoid 
this, it is sufficient to replace the term 
Eq.® by any convenient constant. Then, there is no 
harm in terminating the chain even though the last im- 
age is neither fixed, nor a minimum of the free energy. 
This gives a completely unbiased approach to the deter- 
mination of the MFEP. 



C. Relation to previous approaches 

As far as I am aware, the present use of the NEB 
method for exploring nucleation in the context of DFT 
is novel. However, the problem has been studied using 
other approaches. The properties of the critical cluster 
are in principle accessible since it is an extremum of the 
free energy and several authors have calculated its prop- 
erties, albeit using less sophisticated DFT models than 
that employed here ffi [Til. l27l|. Since the critical cluster is 
a saddle point, rather than a minimum, in the free energy 
surface, it is still difficult to isolate. One would like to 
make an initial guess as to its structure and then to solve 
the Euler-Lagrange equation which inevitably involves an 
iterative refinement of the initial guess until stability is 
reached. However, because the critical cluster is a saddle 
point, and because the only stable minima are the uni- 
form liquid and vapor, one finds that an initial guess that 



is too small evolves towards smaller and smaller clusters 
until the uniform state is reached, while one that is too 
large evolves towards larger and larger clusters until the 
uniform state of the other phase is reached. Oxtoby and 
Evans therefore used this behavior to bracket the critical 
cluster with successive initial guesses and extracted the 
properties of the critical cluster at an intermediate point 
in the iterative calculations [l^- This technique gives cor- 
rect results but is obviously messy and the NEB with the 
climbing image can be seen as an improvement since it 
gives the critical cluster in one run. 

The properties of non-critical clusters have been stud- 
ied by Lee, Telo da Gama and Gubbins[28j and by Talan- 
quer and Oxtoby Q. In the study of Lee et al., finite-sized 
volumes where used. By varying the size of the volume 
and/or the density of material within the volume, clus- 
ters of different sizes become stable - presumably because 
the system wants to phase-separate. The problem in this 
case is that it is not clear what the analogies can be 
drawn between these clusters and those forming in open 
systems. Talanquer and Oxtoby used a method insp ired 
by this approach and developed by Reiss et alpl. 129. [30|. 
In the language used in this paper, the Talanquer-Oxtoby 
method consists of minimizing the grand potential, Cl[n], 
subject to the constraint that the total number of atoms 
in a volume v be fixed at a given number, i, 



p(r)dr. 



(12) 



(The relation between this formulation and that of ref. @ 
is given in the Appendix.) Talanquer and Oxtoby fur- 
thermore adjust the density outside the volume so that 
the density profile is continuous. The idea here is, in 
some sense, to reduce the number of degrees of freedom 
required to explore density function space to just two, 
the parameters i and v, thus reducing the difficulty of 
the problem of finding the MFEP. However, as shown 
in the Appendix, this amounts to changing the chem- 
ical potential so as to stabilize a cluster with the pre- 
scribed number of atoms. In other words, the clusters so 
obtained are just the critical clusters at different chemi- 
cal potentials. These would not appear to be physically 
equivalent to clusters of different sizes at fixed chemical 
potentials. Uline and Corti have used a similar method to 
explore bubble nucleation but without the adjustment of 
the density outside the clusters This circumvents the 
objection raised above, but the method has other prob- 
lems: the generation of discontinuous density profiles and 
the appearance of spurious instabilities 17J . Thus, the 
constraint method, while physically appealing, is not a 
reliable approach for the exploration of density-function 
space. 
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D. Parametrized profiles : Classical nucleation 
theory 

A somewhat simpler alternative to minimizing the 
free energy with respect to a density function is to use 
parametrized density profiles of the form n(r; T) where T 
represents one or more parameters. An example is the hy- 
perbolic tangent profile, n(r; A, zq) — pi + \{p v — Pi)(l + 
tanh(— A(z — Zq))) used to described planer interfaces. 
For a spherical geometry, a closely related form, 



P(r) =Pi+ (pv - Pi) 



1 + br 1 - tanh (A(r - R)) 
1 + (jgrf 1 - tanh (-AR) 



(13) 

has the advantage of having zero derivative at the origin 
and of showing the expected behavior n(r) ~ exp(—ar)/r 
at large r[3l|. This profile will be used as an initial guess 
for the NEB calculation. 

An even simpler parametrization is of particular im- 
portance. Using n(r,R) = pi&{R — r) + p v Q(r — R), 
where the step function Q(x) is one for x > and zero 
otherwise, in conjunction with a simple van der Waals 
free energy functional gives 

Air 

A0SI(R; n, T) = —R (f3P(p v ;T) - 0P(pf,T)) (14) 



4vri? 2 7 c 



Pi - Pv 



Pi 



coex Pv,coex 



where P(p) is the bulk pressure at density p, "f c ,coex is the 
surface tension at coexistence and pi tCoe x and @y <C oex are 
the liquid and vapor densities at coexistence [l7|. This is 
just a slight generalization of the CNT model for the free 
energy of a bubble. (The classical model results from set- 
ting pi,coex-Pv,coex = Pi~Pv as will be done henceforth.) 
This expression for the excess free energy has minima at 
R = and R — > oo corresponding to the uniform vapor 
and liquid, respectively and a maximum at the critical 
radius, 



R c 



27c 



0P(PV,T) - pP(pi;T) 



and 



A0Sl{Rc; fi,T) 



16tt 



7 3 

Icoe.c 



3 (0P( Pv ;T) - 0P{pi;T)f 



(15) 



(16) 



E. Computational details 

In the calculations presented below, the radial variable 
was discretized into 40 points per hard-sphere diameter. 
The profiles were discretized out to a maximum radius of 
20 hard-sphere diameters giving 800 points in total. The 
contributions to the free energy coming from the density 
profiles at larger distances were taken into account by 
assuming the density to be constant and equal to the 
appropriate uniform (liquid or vapor) density at larger 



distances. As alluded to above, the results were very 
insensitive to the value of the spring constant, which was, 
in most cases, fixed at 1 in LJ units. For calculations very 
close to the spinodal, where the energy barriers are very 
small, a smaller spring constant sometimes proved to give 
better convergence. 

The modified Euler-Lagrange equation was solved us- 
ing the fast inertial relaxation engine 32]. In this algo- 
rithm, a fictitious time variable is introduced and each 
component of each image is treated as a dynamical vari- 
able moving in response to the "forces" . The algo- 
rithm involves a quench so that the system relaxes to- 
ward a state of zero force. This algorithm appears to 
be one of the most efficient ways to implement the NEB 
method [33l|. In the present case, it was found to be nec- 
essary to treat the images as being dynamically inde- 
pendent: that is to say, each image had its own time 
variable, time step and quenching variables. The pa- 
rameters governing the quenches were the same as those 
given in ref.[32| and the quench was deemed to be con- 
verged when the root-mean-squared force was less than 
5 x 10~ 5 in LJ units when only the energies were of in- 
terest. The convergence of the profiles, particularly near 
the core, r = 0, are slow because of the r 2 weighting of 
the forces so the algorithm was run until the rms force 
was less than 1 x 10~ 5 when the profiles were desired. 



III. RESULTS 

A. General aspects droplet nucleation 

In this section, the nucleation of liquid droplets from 
a supersaturated vapor will be examined. Before turning 
to specific examples which will mainly involve a com- 
parison to existing simulation data, it is interesting to 
consider some general properties of the present method. 
Two important properties of the critical cluster are its 
size and excess free energy. The size of a cluster may 
be unambiguously defined as the number of atoms in the 
system relative to the number in the metastable state, 



An = (p(r) - Poo)dr 



(17) 



where in the case of droplets (bubbles) poo is the density 
of the bulk vapor (liquid) at the specified chemical poten- 
tial and temperature. Figure Q] shows the size and free en- 
ergy barrier of the critical droplets calculated for the LJ 
potential with cutoff r* = 4 at temperature T = 0.8 and 
excess chemical potential p — p C oex = 0A5p coex , corre- 
sponding to a supersaturation of Sp = P v j ' P coe x = 2.27. 
(Note that with this cutoff, the critical point is calculated 
to occur at T c = 1.25.) The calculations was performed 
using 20 images which were initialized as modified hyper- 
bolic tangents, Eg . (fl"3"|) , with increasing radii. The figure 
shows the results of calculations in which the initial guess 
of the width of the interfaces was very narrow {A = 5 in 
Eq. (|13|) ) which is far from the optimized-sigmoidal value 
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(A ~ 1) as well as the result using initial guesses near 
the optimized-sigmoidal value. It is apparent that, even 
though the initial paths are very different, the final re- 
sults are indistinguishable in shape. Note that one of the 
final curves is longer than another because the choice of 
which image is the "climbing image" , that locates the 
critical droplet, was based on the initial energies and 
these differed in the two runs. Since the climbing im- 
age always ends up at the maximum of the free energy, 
the number of images before and after the maximum can 
differ if the initial curves are sufficiently different. This 
figure also shows that the final free energy barrier is indis- 
tinguishable from the initial optimized-sigmoidal profile 
for large clusters which just confirms that large clusters 
are indeed sigmoidal in shape. However, as shown in the 
inset, the converged barrier is definitely lower than the 
optimized sigmoidal barrier. 

Since the calculations sample the entire nucleation 
pathway, a picture of the development of droplets is 
given. Figure [5] shows the sequence of droplets along the 
MFEP for the same conditions as above. The figure indi- 
cates that the droplets form via a gradual increase in the 
density near the core. This is in contrast to the picture in 
CNT where small droplets are assumed to have the den- 
sity of the bulk liquid over a small volume. In fact, as 
the figure shows, even very small droplets appear to have 
finite volume. This can be quantified by calculating the 
cquimolar radius, defined as the radius of a sphere with 
uniform density p(0) that has the same excess number of 
atoms as the actual density profile, 

Air 

—Rl (p(0) - poo) = An. (18) 

This is shown in Fig. [3] and the results clearly support 
the contention that the radius of the droplets does not 
go to zero. 

B. Comparison to simulation 

There have been many simulation studies of the forma- 
tion of liquid clusters in a supersaturated vapor. Simu- 
lations of stable clusters are possible using constant par- 
ticle number, volume and temperature (NVT): see e.g. 
[H, H3) HE Hfl- Beginning with an unstable density, 
the system will spontaneously phase separate into one 
or more liquid droplets surrounded by vapor. However, 
this does not correspond to nucleation of droplets in the 
laboratory which normally occurs as constant pressure or 
at constant volume but with volumes so large that the 
vapor pressure remains nearly constant. In contrast, the 
NVT simulations have mostly been performed on small 
systems so as to result in the formation of a single droplet 
and the size of the droplet and density of the surrounding 
vapor is then a function of the size of the simulation cell. 
One exception is the work of Oh and Zeng[(| H(| which 
was done using relatively large systems. Comparison to 
their work will be made below. 



Perhaps the cleanest simulations, from the standpoint 
of a comparison to theory, have been those of ten Wolde 
and Frenkel(tWF) carried out in the constant number, 
pressure and temperature - or NPT -ensemble (37) • By 
means of umbrella sampling, they were able to deter- 
mine the properties of unstable clusters of all sizes, at 
least for one value of the supersaturation. Since their 
method gives a faithful sampling of the NPT ensemble, 
direct comparison of theory to simulation is relatively 
straight forward, with some caveats. For example, the 
DFT calculations are performed in the grand ensemble: 
physically, these are expected to give the same results 
as in both cases, droplets are nucleated within a back- 
ground vapor of essentially constant density. More for- 
mally, it is easy to see that the free energy differences 
should be the same in both ensembles 8]. tWF used a 
truncated and shifted potential with cutoff r c = 2.5a and 
all simulations were performed at a reduced temperature 
of T* = k B T/e = 0.741. In the following, the reduced 
density is p* = pa 3 and the reduced length is r* = r/a. 

As stated above, the MC-VDW model was solved using 
the empirical JZG equation of state. However, the JZG 
EOS was developed for an infinite ranged potential. The 
usual means to take account of a finite-ranged potential 
is to introduce mean-field corrections [26] . For relatively 
large cutoffs, this gives an accurate account of the ther- 
modynamics but for cutoffs as short as that used by tWF, 
the corrected equation of state is not very accurate. In 
particular, it gives a critical temperature of T* = 1.0366 
whereas simulation gives T* = 1.085. This difference of 
about 4% is important as the surface tension goes to zero 
at the critical temperature and is thus very sensitive to 
its value. Fortunately, it has been shown that the effect of 
this inaccuracy can be almost completely eliminated by 
invoking the law of corresponding states and comparing 
theory and simulation at equal values of T/T c [38]. This 
is the strategy used here so that the theory is evaluated 
at T/T c = 0.741/1.085 = 0.683 or, given the theoretical 
value of the critical temperature, T* — 0.708. 

Figure [5] shows the excess number of atoms in the crit- 
ical cluster and the free energy barrier as functions of the 
supersaturation, Sp — P v / ' P CO ex where P v is the pressure 
in the bulk vapor and P CO ex is the pressure at coexis- 
tence. The agreement is very good even at quite high 
supersaturations and correspondingly small clusters and 
is seen to be far better than the predictions of CNT. The 
agreement for small clusters is particularly important as 
the clusters are so small, with a radius less than 2a, that 
virtually all atoms in the cluster are part of the interface: 
in other words, the system is extremely inhomogeneous. 
Figure [5] shows the theoretical prediction for the struc- 
ture of the critical cluster at Sp = 1.535 compared to 
the simulation results. In both cases, the structure near 
r = is poorly resolved: in the simulations because there 
are few atoms and poor statistics, and in the calculations 
because the r 2 weighting in the integrals means that this 
region has very little effect on the free energy. Thus, 
ten Wolde and Frenkel only report the density profile for 
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FIG. 1: (Color on line) Initial and final free energy barriers for two different initial conditions: A = 5 and optimized hypertan- 
gent. 




FIG. 2: (Color on line) Profiles of subcritical and supercritical droplets for T* = 0.8 and Sp = 2.27. The system starts in a 
uniform state with density p*(r) = p(r)a 3 = 0.0207 corresponding to an excess number of atoms An = 0. The various curves 
represent various points along the MFEP. The first few curves are labeled with the value of An and the critical cluster is 
indicated with an arrow. 
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FIG. 3: The equimolar radius as a function of the excess number of atoms in a droplet for T* — 0.8 and Sp = 2.27. The 
numbers give An for the indicated clusters. The cluster marked with the arrow is the critical cluster. 
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FIG. 4: (Color on line) Comparison of theory (circles) and simulation results of ten Wolde and Frenkel[37l (squares) for (a) the 
excess number of atoms in the critical cluster and (b) the free energy barrier height as functions of the supersaturation. 
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FIG. 5: (Color on line) Comparison of theory (lines) and the simulation results of ten Wolde and Frenkel[37| (symbols) for the 
structure of the critical cluster for two different values of the supersaturation, S p . 



r* > 1 in ref.[37( although more data is shown in the 
figure here. The theory is clearly very accurate in giving 
the correct overall size and shape of the density profile. 
It does show some structure near the core that is absent 
from the simulation but this is because the theory does 
not take into account the smearing of the interface ex- 
pected to result from both center of mass motion Q and 
capillary waves [3Tjj. 

All of the properties compared so far are restricted 
to the critical nucleus. However, one of the most im- 
pressive aspects of the work of tWF is that they were 
able to determine the free energy barrier as a function 
of cluster size for one particular value of supersaturation 
(S = 1.535). This allows for a check of the novel ap- 
proach used here to determine the MFEP as well as its 
relevance to the nucleation problem. There is, however, 
one subtlety in making this comparison. tWF character- 
ized the barrier as a function of cluster size rather than 
the excess number of atoms in a cluster. Their definition 
of whether or not an atom was in the liquid or vapor 
was based on the local density: an atom was classified 
as part of a liquid cluster if it had at least 4 neighbors 
within a distance of q = 1.5cr. There is no practical, ex- 
act way to translate this into a criteria that can be eval- 
uated theoretically so the following heuristic procedure 
was used. The number of neighbors of within a distance 
q of an atom in the uniform bulk system at density p is 
n = 4-7T Jq pg(r; p)r 2 dr where g(r, p) is the pair distribu- 
tion function(PDF). I have used this expression, together 
with the usual first order Weeks-Chandler- Anderson per- 



turbative approximation for the PDF [13, SO, 0, E2 , to 
determine that n > 4 occurs for p* > 0.32. Hence, in the 
theory, all regions with density satisfying this inequality 
were classified as liquid. For the critical cluster, where 
theory gives and excess number of atoms of 315 com- 
pared to 330 reported by tWF, this procedure gives a 
theoretical cluster size of 265 compared to 285 found in 
the simulations. It would appear that this is a sensible 
way to calculate "cluster size" in the theory. 

Figure [6] shows a comparison of the nucleation barrier 
as a function of cluster size as determined from simula- 
tion and theory. The two are in remarkable agreement 
with the theoretical values about 1.5fcsT smaller than 
those observed in the simulation at the barrier's maxi- 
mum. This is consistent with the fact that the theory 
determines the MFEP whereas the simulations report a 
thermal average which will also involve nearby, higher 
energy, states. This agreement gives strong empirical 
validation of the present theoretical approach. 

Oh and Zeng have performed a set of Monte Carlo sim- 
ulations of large systems at high supersaturation where 
the object was to generate an equilibrium distribution 
of clusters from which the free energy barrier can be 
extracted^ H0|. Their work is complementary to that 
of ten Wolde and Frcnkcl in that it focuses on higher su- 
persaturations, and consequently much smaller clusters. 
However, direct comparison of Oh and Zeng's results to 
theory is complicated by their definition of the supersatu- 
ration which involves a sum over the entire population of 
clusters. In one case, however, they do report the pres- 
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FIG. 6: (Color on line) The nucleation barrier at supersaturation Sp = 1.535 as a function of cluster size. The squares are 
from simulation of ten Wolde and Frenkel[13| and the circles are from the theory. 
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sure of the ambient vapor as Pa 3 / (fcsT) = 0.01166 at 
T* = 0.67[36|. Since the potential cutoff used in the 
simulations was r c = 4.5<r (with no shift), it can be 
assumed that the JZG equation of state with a mean- 
field correction for the cutoff is essentially exact and no 
corresponding-states adjustments are necessary. In this 
case, the main uncertainty comes from the definition of 
the cluster size, denoted no- Oh and Zeng define two 
atoms to be in the same cluster if they are within 1.5(7 
of one another. The theoretical cluster size was therefore 
calculated as in the case of tWF except the limiting den- 
sity was chosen to be that at which there is less than one 
neighbor in a sphere of radius 1.5cr. Figure [7] shows the 
free energy barrier determined from simulation compared 
to that calculated here. In particular, the barrier height 
is calculated to be 6.5ksT compared to the value from 
simulation of 8.2ksT while the peak occurs at almost the 
same place, no ~ 20. Although the relative error in the 
barrier height is greater than in the case of the tWF data, 
the comparison is actually remarkably similar in that the 
theoretical barrier is lower than the observed barrier by 
almost exactly the same amount found in the compar- 
ison to tWF while the cluster size is virtually identical 
to that found in the simulations, again as found in the 
comparison to tWF. Given the very small cluster size, 
the uncertainty in the thermodynamic state and in the 
estimate of the cluster size, and the fact that the simula- 
tions involve many clusters while the calculations are for 
an isolated cluster, the agreement is probably as good as 
could be expected. 

The comparison to the results of Zhukovitskii [43| is 
more problematic for two reasons. First, there was no 
fixed cutoff in the simulations: the potential was essen- 
tially infinite-ranged but the volume is finite so that the 
number of neighbors a given atom interacted with de- 
pended on how close that atom was to the (spherical) 
cell boundary. The radius of the cell was Rz = 8cr so the 
effective cutoff is in the range < r c < 8a. Zhukovitskii 
reports the density of the vapor at coexistence as being 
Pv,coex = 0.00243 so I have adjusted the cutoff to approx- 
imately match this, which happens at r c ~ 5<7. Clearly, 
this is only an approximation and probably represents the 
greatest source of error. Next, Zhukovitskii reports the 
supersaturation in terms of the variable S p = p v /p VyC oex 
however, this was only estimated assuming the vapor is 
an ideal gas: apparently, the actual vapor pressure was 
typically about 4% lower than this[43j], so in the calcula- 
tions, it was assumed that p v — 0.95S p p VtC oex- With such 
a large cutoff, the error in the equation of state should 
be minimal so there was no need for the kinds of adjust- 
ment required for the tWF cutoff. Given these caveats, 
Fig. [8] shows the size of the critical cluster as a function 
of supersaturation as reported by Zhukovitskii compared 
to the calculations. Given the uncertainty in the cutoff 
and in the actual density of the vapor, not to mention the 
fact that the simulation does not correspond precisely to 
any standard ensemble, the agreement is satisfactory. 



C. Nucleation of Bubbles 

The same methods can be used to study the nucleation 
of bubbles in a superheated liquid. Figure [9] shows the 
height of the nucleation barrier as a function of bubble 
size for different supersaturations for T* — 0.8 and a 
cutoff of r* = 4.0. Note that for bubbles, An is negative 
and the supersaturation is given here as S u = A1 ~ Alc ° ex . 
. In contrast to recent claims that there is an activated 
instability in bubble nucleation fl~oj] . the results indicate 
continuous paths for a wide range of supersaturations. 
It is possible that at high absolute value of supersatura- 
tion the free energy is concave near An = but there 
is no sign of a non-classical instability, even when the 
free energy barrier is less than 2fcsT. Figure [TU] shows 
a sequence of bubbles for S p = —0.15. As in the case 
of droplet nucleation, the width is apparently always 
greater than zero with the bubble nucleating via an ini- 
tial, gradual lowering of the central density followed by 
a slow broadening into a typical sigmoidal shape. Note 
that the density at the center of the critical cluster shown 
in Fig. [TU]is more than half the liquid density and, hence, 
much greater than the density of the gas phase being nu- 
cleated. This represents a large difference from the as- 
sumption of CNT that the critical cluster has the density 
of the bulk vapor inside the bubble. 

IV. DISCUSSION 

The nucleation of a stable phase from a metastable 
phase is best understood as a transition between two local 
minima in the free energy surface. As such, it is concep- 
tually the same as any problem that involves the crossing 
of a free energy barrier and, in particular, bears strong 
similarities to the description of chemical reactions. It 
is therefore not surprising that the methods used to de- 
termine reaction pathways can be usefully applied to the 
problem of nucleation. 

There are in fact many methods used to study reac- 
tion pathways including eigenvalue-following [II H], the 
string method 14511 . the determination of the maximum 
likelihood path[46[ among others, including the NEB 
method used here. The NEB method may not be the 
optimal method for nucleation, but it has the advantage 
of being very simple to implement and of being robust. 

There has been much discussion recently concerning 
the possibility of non-classical mechanisms of liquid- 
vapor and of vapor-liquid nucleation. In particular, 
Bhimalapuram, Chakrabarty and Bagchi have observed 
a spinodal-like breakdown at high supersaturations in 
Ising model simulations of condensation[47| although 
the significance of their observation has recently been 
questioned [48j . Similarly, Uline and Corti have claimed 
a similar type of behavior for boiling based on a combina- 
tion of simulations and DFT calculations 10]. A similar 
observation concerning crystallization has been made by 
Parinello 49] . What all of these studies have in common 
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FIG. 9: (Color on line) The free energy barrier for bubble nucleation as a function of bubble size at several different values 
of supersaturation as calculated from the theory. The curves are labeled with their supersaturation, The inset shows an 
expansion of the small An region. 
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FIG. 10: (Color on line) A sequence of bubble profiles along the MFEP for supersaturation S M = —0.20. The size, in terms of 
atomic deficit or —An, is given for several profiles and the critical profile is marked with an asterisk. 
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is the observation of nonstandard elements in the free 
energy surface: either an additional local minimum (47j 
or some type of discontinuity [l^. The present work sup- 
ports none of these observations. In all cases, even at 
very high supersaturation, the free energy surface ap- 
pears to be well-defined and to possess no non-classical 
minima. For the nucleation of droplets from a vapor, 
these results are consistent with the simulations of ten 
Wolde and Frenkel[37| and of Oh and ZengQ, neither of 
which saw evidence of non-classical features in the free 
energy barrier. In the case of the nucleation of bubbles in 
a superheated liquid, the calculations may indicate that 
at very high supersaturations, such that the free energy 
barrier is less than 2ksT, the free energy barrier is not a 
convex function of its size, but it is not clear that this is 
of any significance. Furthermore, there is evidence, based 
on simple analytically tractable models, that the discon- 
tinuities observed by Uline and Corti are an artifact of 
the constraint method used to explore the free energy 
surface jl7j|. Thus, while the present results cannot de- 
cisively settle the question one way or another, they do 
tend to support the classical picture of liquid-vapor nu- 
cleation as a simple matter of barrier crossing. 
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APPENDIX A: THE CONSTRAINT METHOD 

The constraint method described by Talanquer and 
Oxtoby consists of demanding that the number of atoms 
within the volume v be fixed at some value i. Further- 
more, the density outside the volume is prescribed to be 
a fixed value, poo. Minimizing the grand potential un- 
der these constraints can be formulated using Lagrange 
multipliers so that one minimizes the functional 



F[n] - (J, J n(r)dr - a ( / n(r)dr - i ) (Al) 

'r<R 



7(r) (n(r) - p^) dr 



r>R 



with respect to n(r), a, and 7(r). For r < R this gives 



SF[n] 
5n 



p — a 



(A2) 



0=1 n(r)dr — i 

'r<R 



so that it is clear that the effect of the Lagrange multiplier 
is simply to shift the chemical potential. The resulting 
density profile is just that of the system at the shifted 
chemical potential. To make contact with the work of 
Talanquer and Oxtoby, it is necessary to separate the 
functional derivative occurring in this expression into its 
ideal gas and excess contributions by writing 



8F[n] 
Sn 



fc B Tlogn(r) + p ex {r) 



(A3) 



This allows the first of the variational equations to be 
rearranged to give 



n(r) = exp (/3p + f3a - j3p ex {r)) 



(A4) 



Substitution into the second variational equation gives 
an expression for the Lagrange Multiplier 

= exp(/3a) / exp (fip - (3p ex {v)) dr - i 

Jr<R 



(3a + f3p = log(z) -log exp(-(3p ex (r))drj (A5) 

Thus, the equation for the density, for the case r < R, 
becomes 



60F[r 
Sn 



log(i) - log 



cxp(-f3p ex (r))dr) (A6) 



r<R 



which is, essentially, Eq. 17 of ref. ([§]). 
The variational equations for R < r are 

= — fi- 7 (r) 

on 

= n(r) - poo 



(A7) 



for very large r, far away from the interface, it will be 
the case that 



V 



-i fSF[r 



\ 5n 



n=p 



df(Poc) 
dpoo 



(A8) 



so that in this region, p + j(r) 



where 

the last term means the chemical potential corresponding 
to the density p^. In particular, if R is sufficiently large 
that the cluster interface is located at much smaller val- 
ues of the radius, as appears to be the case in the work of 
Talanquer and Oxtoby, then a continuous density profile 
will basically require that a + p, = 7(r) + ji = p(poo) so 
that the procedure simple amounts to an overall shift of 
the chemical potential. 
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